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Abstract. 

This  paper  introduces  a  computational  method  for  use  with  eddy  current  damage  detection  tech¬ 
niques.  To  identify  the  geometry  of  a  subsurface  damage,  an  optimization  algorithm  is  employed 
which  requires  solving  the  forward  problem  numerous  times.  In  order  for  the  method  to  be  effective 
in  a  practical  setting,  i.e.,  in  real-time  applications,  the  forward  algorithm  must  be  extremely  fast  and 
accurate.  Therefore,  we  have  chosen  an  approach  based  on  reduced  order  Proper  Orthogonal  Decom¬ 
position  (POD)  methods.  This  allows  one  to  create  a  set  of  basis  elements  using  snapshots  with  either 
numerical  simulations  or  experimental  data.  The  data  is  organized  in  an  optimal  way  allowing  one  to 
use  a  reduced  number  of  basis  elements,  resulting  in  a  fast  algorithm  while  still  obtaining  an  accurate 
approximation  to  the  solution.  We  first  derive  the  model  associated  with  the  chosen  nondestructive 
evaluation  (NDE)  technique  and  prove  some  well-posedness  results  for  the  model.  We  then  introduce 
the  proposed  computational  methodology  and  test  it  on  both  numerically  simulated  data  as  well  as 
experimental  data  obtained  from  a  GMR  (Giant  Magnetoresistive)  sensor.  The  results  demonstrate 
that  the  method  is  extremely  efficient  and  accurate. 


1  Introduction  and  Problem  Formulation 

As  technology  continually  advances,  the  field  of  nondestructive  evaluation  is  in  continual  need  of  new 
techniques  and  instruments  to  improve  the  accuracy  and  efficiency  of  locating  and  characterizing  sub¬ 
surface  damages.  We  attempt  to  develop  a  new  methodology  which,  when  coupled  with  already  existing 
techniques,  can  help  decrease  the  total  computational  time  required  to  detect  and  explicitly  characterize  a 
damage  within  a  material.  This  is  necessary  in  practical  settings  where  the  methods  must  be  fast  and  ac¬ 
curate,  producing  real-time  results.  Given  data  obtained  from  a  measuring  device  such  as  the  GMR  (Giant 
Magnetoresistive)  sensor,  we  seek  to  locate  and  parameterize  the  damage  while  minimizing  the  amount  of 
time  required  to  complete  this  task.  To  this  end,  we  formulate  and  develop  an  appropriate  model  used  in 
describing  the  variation  in  the  data  as  a  function  of  a  damage  within  the  sample  and  present  computational 
methods  along  with  numerical  results  to  support  the  efficacy  of  our  approach. 

1.1  Description  of  the  Test  Problem 

An  advanced  method  of  damage  detection  uses  a  device  such  as  the  GMR  sensor  in  the  context  of  eddy 
current  methods  [1,  2].  In  a  standard  experimental  setting,  as  depicted  in  Figure  1,  a  thin  conducting  sheet 
carrying  a  uniform  current  is  placed  above  (or  below)  the  sample.  The  current  within  the  sheet  induces  a 
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magnetic  field  perpendicular  to  it  that  in  turn  produces  a  current  within  the  sample,  called  an  eddy  current. 
When  a  flaw  is  present  within  the  sample,  the  flaw  disrupts  the  eddy  current  flow  near  the  flaw  and  this 
disturbance  is  manifested  in  the  magnetic  flux  density  detected  by  the  measuring  device. 


Conducting  Sheet 


Sample  Material  Flaw 


Figure  1:  Inspection  Process  using  a  SQUID 

For  illustrative  purposes,  we  will  assume  uniformity  along  the  width  (z  direction)  of  the  sample, 
therefore  reducing  the  three-dimensional  problem  to  a  two-dimensional  problem.  To  test  the  feasibil¬ 
ity  of  reconstructing  the  geometry  of  the  damage,  we  consider  the  damage  (which  we  shall  refer  to  as 
a  “crack”)  to  be  rectangular  in  shape.  In  other  words,  we  assume  the  crack,  located  at  a  certain  depth 
within  the  sample,  has  a  fixed  length  and  thickness.  To  further  simplify  the  test  problem,  we  disregard 
the  boundary  effects  of  the  materials  in  the  x  direction  (sample  length)  by  assuming  an  infinite  sam¬ 
ple  and  conducting  sheet  in  that  direction.  Because  we  are  considering  materials  of  infinite  extent,  we 
will  construct  our  forward  problem  by  focusing  on  a  small  “window”,  called  our  computational  domain 
Q  =  {(x,  ?/,  2;)  G  :  0mm  <  x  <  50mm,  —35mm  <  y  <  35mm},  centered  such  that  the  left  boundary 
of  the  computational  domain,  at  location  x  =  0,  is  positioned  in  the  center  of  the  crack  in  the  x  direction, 
i.e.,  the  crack  is  symmetric  through  the  yz  plane  at  x  =  0.  A  schematic  of  the  resulting  two-dimensional 
problem  is  depicted  in  Figure  2  where  it  is  assumed  that  the  sample  (which  is  20mm  thick)  is  composed 
of  aluminum,  the  conducting  sheet  (which  is  0.1mm  thick)  is  made  up  of  copper  and  the  crack  is  centered 
in  the  y  direction  around  the  center  of  the  sample  (i.e.,  around  y  =  —10mm). 


y, 


Air 

/  Conducting  Sheet  (Copper) 

— 1  Sample  (Aluminum) 

50mm 

Air 


Figure  2:  2-D  Schematic  of  Problem 
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1.2  Resulting  Equations  for  the  Test  Problem 


In  our  computational  efforts,  we  employ  the  use  of  the  software  package  Ansoft  Maxwell  2D  Field  Sim¬ 
ulator.  Therefore,  our  equations  are  formulated  to  correspond  to  those  used  by  the  software.  (For  a  full 
derivation  of  the  resulting  equations,  see  [3].)  Although  the  GMR  sensor  measures  the  magnetic  flux  den¬ 
sity  above  the  sample  at  a  certain  location,  our  equations  are  formulated  in  terms  of  the  magnetic  vector 
potential  A,  where  all  the  field  quantities  are  assumed  to  be  phasor  quantities  [4,  5,  6].  However,  one  can 
easily  obtain  the  resulting  magnetic  flux  density  B  through  the  relationship  B  =  V  x  A. 

Using  Maxwell’s  equations  in  conjunction  with  Ohm’s  law  and  constitutive  laws,  we  obtain  an  equation 
for  the  magnetic  vector  potential  A  given  by 


V  X 


1 


V  X  A(x,  y)  =  {a{x,  y)  -f  iuje){A{x,  y)  -  V<?i)  Vx,  ?/  €  G, 


(1) 


where  g  represents  the  magnetic  permeability,  a  represents  the  conductivity,  u)  is  the  angular  frequency, 
and  (j)isa  scalar  potential. 

Since  Eq.(l)  contains  two  unknowns,  A  and  an  additional  equation  is  needed  to  uniquely  determine 
the  solution  for  A.  For  this  we  use  an  integral  constraint  given  by 


/  Jt  •  nda  =  / 

Jcs  Jcs 


{a{x,  y)  -|-  iix!e{x,  y)){A{x,  y)  —  Vf)  ■  nda 


(2) 


between  the  total  current  Ics  flowing  in  the  conducting  sheet  (cs)  and  the  total  current  density  Jt  within  the 
conducting  sheet.  This  is  the  second  equation  used  in  the  software  package  Ansoft  Maxwell  2D  Field  Sim¬ 
ulator  which  we  use  in  our  computational  efforts.  However,  this  only  gives  us  equations  which  completely 
describe  the  magnetic  vector  potential  and  electric  scalar  potential  in  the  conducting  sheet.  Nonetheless, 
the  conducting  sheet  is  the  only  region  in  which  a  source  current  of  the  form  .Jg  =  —aV(j)  is  present 
and  hence  is  the  remaining  regions  we  intuitively  assume  there  is  no  change  in  potential,  in  other  words 
=  0.  This  gives  us  appropriate  equations  in  which  the  magnetic  vector  potential  A  can  be  uniquely 
determined  if  appropriate  boundary  conditions  on  A  are  specified. 

In  Section  1 . 1  we  assumed  the  sample  was  of  infinite  extent  with  the  crack  being  symmetric  in  the  x 
direction.  In  other  words  on  the  x  boundaries,  we  assume  the  fields  on  both  sides  of  the  boundary  oscillate 
in  the  same  direction.  To  account  for  the  even  symmetry,  we  assign  Neumann  boundary  conditions  to  these 
boundaries.  In  a  similar  manner,  we  assume  the  y  boundaries  are  “sufficiently  far”  away  from  the  sample 
and  scanning  area  so  as  to  not  effect  the  overall  measurements.  Indeed,  as  one  moves  farther  away  from  the 
sample  and  conducting  sheet,  the  magnetic  vector  potential  A  tends  to  zero.  Therefore,  on  the  y  boundaries 
we  assign  null  Dirichlet  boundary  conditions.  The  magnetic  vector  potential  A  is  thus  determined  by 


=  {a{x,  y)  -f  iue{x,  y)){-iuA{x,  y)  -  Vf)  Vx,  ?/  €  D, 

(3) 

If.g  —  Jt  -  nda  — 

J  {a{x,  y)  -\-  iix!e{x,  y)){—iix!A{x,  y)  —  V<?i)  -  nda 

(4) 

and 

=  0  yx,y^fi\cs 

(5) 

with 

A(x,-35)  =  0  =  A(x,35) 

VA  -111(0,2;)  =  0  =  VA  -111(50,2;). 
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2  Well-Posedness 

In  this  section  we  consider  the  existence  and  uniqueness  of  a  weak  solution  to  the  above  boundary  value 
problem  on  a  general  domain  given  by 

^  2/;  '2')  £  •  ^min  ^  ^  ^  ^maxt  Umin  —  U  —  Umax} 

for  which  our  test  problem  is  a  specific  example.  Let  i?  =  L2{D)  and  L  =  {xp  ^  '^{x-,ymin)  =  0  = 

Umax)}  where  we  use  the  standard  Sobolev  space  notation,  =  {xp  e  L?{D)  :  Vxp  e 

and  note  that  we  interpret  pointwise  evaluation  of  functions  (along  the  boundary  and  elsewhere)  in  terms  of 
a  trace  operator  for  which  we  suppress  notation  throughout  this  paper  [7] .  We  denote  by  {(p,ip)  =  (pxpda 
the  standard  inner  product  in  H  and  xp)v  =  Jq  Vp)  ■  V'tpda  the  (77^ -equivalent)  inner  product  in  V. 

We  note  that  in  this  two-dimensional  problem,  the  term  Vp)  can  be  proven  to  be  piecewise  constant  as 
done  in  [3].  In  doing  this,  Vp)  can  be  written  in  terms  of  the  magnetic  vector  potential  in  the  conducting 
sheet  by  solving  for  in  (4).  We  can  then  reduce  the  system  into  an  integro-differential  equation.  Using 
integration  by  parts  together  with  natural  boundary  conditions  and  imposed  conditions  on  test  functions 
xp  ^V,  the  variational  form  is  given  by 

( Vvl,  S/xp)  +  {PiA,  xp)  +  P2  /  Ada  /  xpda  =  /  fxpda.  (6) 

J  cs  J  cs  J  cs 

where  Pi  =  icoir{a  +  icoe),  P2  =  -  ^  and  /  = 

We  consider  both  the  existence  and  uniqueness  of  the  solution  A  to  (6)  as  well  as  the  continuous 
dependence  on  the  parameters  which  represent  the  damage  in  the  context  of  a  Gelfand  triple  setting  V  ^ 
77  ~  77*  ^  U*  where  we  have  that  the  embedding  U  ^  77  is  dense  and  continuous  with 

\'4’\h  <  k\'4^\v  for  all (7) 

where  the  norm  in  V  will  be  denoted  by  |  •  |y  and  |  •  |  will  denote  the  norm  in  77  for  the  rest  of  this  section. 

Using  this  notation,  we  have  the  following  theorem.  For  full  details  and  proof,  we  refer  the  reader  to 
[3,  8]. 

Theorem  2.1.  There  exists  IF  =  IF  {/i,  e,  U)  such  that  for  source  frequencies  fs  <  lF{ii,  e,  U),  there  exists 
a  unique  weak  solution  A  to  (6). 

Furthermore,  we  let  q  represent  the  four  comers  of  any  quadrilateral  damage,  i.e., 
q  =  [(xi,j/i),  (2:2, 2/2),  (2^3, 2/3):  {^4,2/4)]  is  a  vector  in  x  R^  or  equivalently  R*.  We  denote  by  Qad 
the  set  of  admissible  parameters  q  where  it  is  assumed  Qad  is  a  compact  subset  of  R® .  Then  for  any  two 
damages  given  by  parameters  q  and  q  in  Qad,  let 

=  l!q-  q||  =  [{^I  -  xff  +  {yi  -  yff  +  ...  -f  {xi  -  x^f  +  {yi  -  yif]^'‘^  (8) 

to  be  the  standard  Euclidean  norm  in  R* .  We  denote  by  5D  (see  Figure  3)  the  points  in  U  which  are  either 
in  the  damage  represented  by  q  or  q  but  which  are  not  in  both.  In  other  words,  let  Uq  represent  the  points 
{x,  y)  in  U  within  the  damage  given  by  q  and  Uq  the  points  (x,  y)  in  U  within  the  damage  given  by  q. 
Then  =  Uq  U  Uq  —  Uq  fl  Qq. 

Using  the  terminology  above,  we  have  the  following  theorem.  Again,  we  refer  the  reader  to  [3,  8]  for 
details  and  proof. 

Theorem  2.2.  Assume  the  admissible  parameter  set  Qad  A  a  compact  subset  o/R*.  Then  there  exists 
T  =  IF  {p,  e,  U)  such  that  for  source  frequencies  fs  <  lF{p,  e,  U),  q  — >■  A(q)  is  continuous  from  Qad  to  V. 
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Conducting  Sheet 


damage  q 


damage  q 


Sample 


Figure  3:  The  Area  Represented  by  6Ct 


3  Computational  Method 

To  enable  the  teehniques  used  in  nondestruetive  evaluation  to  be  implemented  in  a  praetical  setting,  we 
must  not  only  locate  and  characterize  subsurface  damages;  we  need  to  do  so  in  a  fast  and  efficient  manner. 
To  develop  a  fast  and  efficient  forward  algorithm,  we  use  the  reduced  order  Karhunen-Loeve  or  Proper 
Orthogonal  Decomposition  (POD)  methodology.  A  unique  feature  of  the  POD  technique  is  its  ability 
to  create  an  ordered  basis  that  models  experimental  or  simulated  data,  capturing  most  of  the  important 
aspects  of  the  data  in  the  first  few  elements.  Typically  the  number  of  basis  elements  required  in  the  forward 
algorithm  is  reduced,  resulting  in  the  development  of  a  faster  algorithm  that  still  maintains  the  accuracy  of 
traditional  finite  element  algorithms. 

3. 1  The  POD  Method 

In  this  section  we  discuss  the  concepts  and  method  of  implementation  of  the  POD  method  in  the  context 
of  damage  detection.  For  details  on  the  general  POD  method,  we  refer  the  reader  to  [9,  10,  11,  12,  13,  14, 
15,  16,  17,  18,  19,  20]  and  the  extensive  list  of  references  contained  therein.  In  the  process  of  detecting  a 
subsurface  damage,  a  device,  such  as  the  GMR  sensor,  is  scanned  above  (or  below)  a  sample  and  data  is 
taken.  The  variation  in  the  field  due  to  a  damage  is  manifested  in  the  data  taken.  Therefore,  in  forming  a 
reduced  basis,  we  want  to  incorporate  the  effects  of  a  damage  in  the  reduced  basis.  This  is  accomplished 
by  taking  “snapshots”  of  the  data  across  various  damages.  In  other  words,  let  q  be  a  vector  parameter 
characterizing  physical  properties  of  the  damage  such  as  length,  thickness,  depth,  etc.  of  the  damage.  Then 
given  an  ensemble  of  damages  vve  obtain  corresponding  solutions,  {A(qj)}^i,  of  the  boundary 

value  problem,  for  magnetic  vector  potentials  which  we  call  our  “snapshots”.  We  then  use  these  snapshots 
to  generate  a  basis  incorporating  the  properties  of  the  various  damages.  (Without  loss  of  generality,  we 
will  denote  the  vector  A  by  its  scalar  nonzero  component  A,  i.e.,  the  component  of  A.) 

The  formation  of  the  POD  basis  can  be  summarized  in  a  few  steps.  In  order  to  successfully  reduce  the 
number  of  basis  elements  while  maintaining  accuracy,  a  single  basis  element  must  contain  aspects  of  each 
damage  q^.  Therefore,  as  explained  in  [20,  21],  we  seek  basis  elements  of  the  form 

Ns 

j=i 

where  the  coefficients  Vi{j)  are  chosen  such  that  each  POD  basis  element  i  =  1, 2, ...,  Ng,  maximizes 
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subject  to  ($j,  =  1-  Forming  the  basis  elements  in  such  a  way  assures  that  a  single 

basis  element  will  contain  information  from  each  of  the  snapshots.  From  standard  arguments,  it  can  be 
seen  that  the  coefficients  Vi{j)  can  be  found  by  solving  the  eigenvalue  problem  CV  =  \V  where  C  is 
given  by 

[C]ij  =  (11) 

Hence,  the  next  step  in  forming  the  POD  basis  is  to  form  the  covariance  matrix  C  and  find  the  asso¬ 
ciated  eigenvalues  and  eigenvectors.  Since  C  is  Hermitian  positive  semi-definite,  we  know  it  possesses  a 
complete  set  of  orthogonal  eigenvectors  and  corresponding  nonnegative  eigenvalues.  In  forming  the  POD 
basis,  we  want  to  be  able  to  readily  decide  which  basis  elements  to  use  in  the  reduced  basis.  To  this  end, 
we  order  the  eigenvalues  along  with  their  corresponding  eigenvectors  such  that  the  eigenvalues  are  in 
decreasing  order, 

Ai  >  A2  ^  ...  >  AjVs  ^  0.  (12) 

We  then  normalize  the  eigenvectors  corresponding  to  the  rule 


1".  •  1".  = 


% 


(13) 


Consequently,  the  POD  basis  element  is  defined  by  Eq.  (9)  where  Vi{j)  represents  the  component 
of  the  eigenvector  of  C. 

We  now  have  the  full  POD  basis  and  need  a  criterion  to  decide  how  many  basis  elements  are  re¬ 
quired  to  accurately  portray  the  data.  In  other  words,  we  want  to  choose  N  such  that  ^  ^ 

span{A{c{j)}^A-  Iii  choosing  this  number  N,  we  compute 


N  Ns 


f=i  f=i 


(14) 


which  represents  the  percentage  of  “energy”  in  span{A{qj)}^A  that  is  captured  in  span{^j}^^^.  Then 
the  reduced  basis  consists  of  the  first  N  POD  basis  elements  where  N  is  chosen  to  capture  the  desired 
amount  of  energy. 

We  now  note  that  span{^i}^A  —  span{A{qj)}^A-  Indeed,  given  any  ^(q^),  we  have 


Ns 

k=l 


(15) 


where 

o^ki^ij)  =  {A{qj),^k)L^{n,'c)  (16) 

as  are  orthonormal  in  L^(0,  (D). 

However,  to  complete  the  analysis,  one  needs  to  be  able  to  calculate  ^^(q)  where  q  is  a  given  param¬ 
eter  not  in  the  set  {q^}^!.  To  this  end,  we  extend  the  approximation  formula  to  obtain 

k=l 


(17) 
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where  Q;A;(q)  can  be  evaluated  through  different  methods.  In  [3],  we  explored  two  different  methods 
one  might  choose:  a  POD/Galerkin  method  or  a  POD/Interpolation  method.  However,  in  this  paper,  we 
only  present  results  in  which  we  use  the  POD/Interpolation  method  with  linear  interpolation  for  the  one- 
parameter  simulated  results  and  cubic  spline  interpolation  for  all  other  results.  For  details  on  these  meth¬ 
ods,  we  refer  the  reader  to  [22]  and  [23]. 

3.2  Simulated  Results 

In  [4],  we  performed  several  trials  in  which  we  assumed  we  had  access  to  various  types  of  data,  such  as 
the  A  field  or  the  B  field  at  various  points  (xj,  yj)  in  O.  We  compared  and  contrasted  the  accuracy  to 
which  we  could  estimate  the  length  I  of  the  damage  based  on  whether  the  A  field  or  B  field  was  used  and 
whether  we  considered  the  field  along  a  single  line,  multiple  lines  or  within  the  entire  region  (which  is  not 
experimentally  possible  and  was  only  tested  for  initial  comparisons).  From  the  results  presented  in  [4],  we 
concluded  that  extremely  accurate  results  were  obtained  only  when  the  y  component  B2  of  the  magnetic 
flux  density  was  used  in  the  cost  criterion,  i.e.,  when  we  used 

^  n  m 

'^(q)  =  9  XI XI 1 10*^^ {Xi,yj,q)  -  10^B2{xi,yj,q*)\'^  (18) 

i=l  j=l 

where  10*  is  a  scaling  factor  accounting  for  the  low  order  of  magnitude  of  the  field  (B2  is  on  the  order 
of  10“*FF()/m),  B^{q)  =  V  x  ^^(q)  is  the  reduced  order  POD  approximation  in  which  ^^(q)  is 
given  by  (17),  and  132  is  “data”  from  a  sample  we  wish  to  characterize.  (In  this  section,  B2  is  obtained 
from  Ansoft  finite  element  simulations  to  which  noise  has  been  added  in  the  usual  manner  (see  [3,  21]). 
Furthermore,  performing  multi-line  scans  or  using  full  region  data  improved  the  results  only  marginally 
and  hence  did  not  warrant  the  extra  effort  and  time  in  collecting  more  extensive  data  sets.  Consequently 
the  results  presented  in  this  section  involve  only  the  least  squares  difference  in  the  B2  field  given  by  (18) 
along  a  single  line  located  1mm  above  the  conducting  sheet. 

In  the  sample  results  presented,  we  focus  on  estimating  a  single  parameter  of  the  crack  assuming 
all  other  parameters  are  fixed  or  estimating  two  parameters  with  one  parameter  fixed.  The  results  below 
summarize  the  feasibility  of  determining  the  length  and  depth  of  a  crack.  In  a  specific  trial  run,  ten  different 
data  sets  (exact  data  with  ten  different  sets  of  added  random  noise)  are  used  where  the  relative  noise  is 
chosen  at  a  10%  noise  level  with  a  confidence  level  of  99.7%  (3  standard  deviations).  (Details  can  be  found 
in  [3]  and  [4].) 

In  determining  the  length  of  the  damage  alone,  we  first  followed  the  steps  outlined  in  Section  3.1  and 
generated  an  ensemble  of  “snapshots”.  Keeping  the  crack  fixed  at  a  depth  of  9mm  with  a  thickness  of 
2mm,  we  varied  the  crack  length  I  from  0mm  to  4mm  in  increments  of  0.2mm.  Using  the  snapshots, 
{A{lj)}^T^‘^^,  we  formed  the  POD  basis  and  determined  that  99.99%  of  the  energy  of  the  system  was  cap¬ 
tured  with  a  single  basis  element  [4].  However,  in  the  specific  trial  runs,  we  found  that  using  4  POD  basis 
elements  improved  the  parameter  estimation  with  data  containing  no  noise,  while  there  was  no  significant 
difference  between  the  use  of  5  basis  elements  over  the  use  of  4  basis  elements.  Therefore,  the  results  are 
based  on  the  use  of  4  basis  elements  in  the  algorithm. 

To  test  the  inverse  methodology,  we  first  try  to  identify  the  length  of  the  damage,  I*  =  1.3mm.  Using 
data  on  a  single  line  above  the  conducting  sheet  with  10%  relative  noise,  an  average  estimated  length  of 
1.2977mm  was  obtained  with  a  variance  of  0.3237  x  10“^mm^  across  the  ten  trials. 


Electromagnetic  Interrogation  Techniques  for  Damage  Detection 


9 


Proceeding  as  we  did  in  estimating  the  length  of  a  damage,  to  estimate  the  depth  of  the  damage, 
we  fixed  the  thickness  of  the  crack  at  0.5mm  with  length  2mm  and  varied  the  depth  of  the  crack  d  from 
0.25mm  to  16.25mm  in  increments  of  0.5mm.  Using  5  basis  elements,  we  were  able  to  accurately  estimate 
a  crack  depth  of  8mm  in  the  presence  of  10%  relative  noise;  an  average  depth  of  8.0631mm  was  estimated 
with  a  variance  0.1180  x  10“^mm^.  Thus,  as  in  the  case  of  estimating  the  length  of  a  crack,  we  can  also 
recapture  the  depth  of  a  crack  quite  accurately  and  efficiently. 

Finally,  we  estimated  both  length  and  depth  simultaneously.  However,  first,  in  [3,  8],  we  discussed  in 
detail  a  need  to  modify  the  assumptions  made  in  the  original  test  problem  to  more  accurately  describe  the 
behavior  of  experimental  data  used  in  the  next  section.  In  short,  the  computational  domain  was  expanded 
beyond  the  edges  of  the  sample  and  snapshots  were  taken  of  the  magnetic  flux  density  data  on  a  single 
line  above  the  conducting  sheet  (instead  of  the  whole  region).  Recall  that  for  the  previous  trials,  we  took 
snapshots  of  the  magnetic  vector  potential  for  the  entire  computational  domain  even  though  in  the  inverse 
problem  we  only  considered  those  data  points  along  a  single  line.  Furthermore,  we  considered  data  across 
the  entire  length  of  the  sample,  instead  of  just  half  the  sample  as  done  previously.  (For  more  details  ,  see 
[3,  Chapter  6].)  As  a  result,  we  also  implement  these  changes  in  the  two-parameter  estimation  problem. 

We  proceed  as  in  the  previous  estimation  problems  by  first  generating  an  ensemble  of  damages.  We 
consider  damages  with  depths  ranging  from  1mm  to  11mm  in  increments  of  2mm  in  combination  with 
lengths  from  0.5cm  to  3.5cm  in  increments  of  1cm  (we  now  consider  longer  damages  similar  to  those  used 
in  obtaining  experimental  data).  We  keep  the  thickness  fixed  at  1mm.  A  total  of  24  snapshots,  {B2{di^  Ij)}, 
f  =  1,  ...,6,  j  =  1,  ...,4  were  generated  using  Ansoft.  Using  4  basis  elements  with  10%  relative  noise 
added,  we  estimated  a  damage  with  depth  d*  =  4mm  and  length  T  =  2cm.  We  obtained  average  es¬ 
timates  of  d  =  4.0855mm  with  variance  0.5516  x  10“^mm^  and  length  I  =  1.980cm  with  variance 
0.9632  X  10“^mm^.  Thus,  even  when  we  estimated  two  parameters  simultaneously,  the  computational 
method  proposed  produced  extremely  accurate  results. 

3.3  Experimental  Results 

Since  the  results  using  simulated  data  were  extremely  promising,  we  designed  an  experiment  in  which 
we  tried  to  detect  and  parameterize  a  damage  within  an  aluminum  sample  using  a  giant  magnetoresistive 
(GMR)  sensor.  The  sample  was  constructed  of  17  layers  of  1mm  thick  aluminum  plates  with  a  slice  cut 
out  of  one  of  the  layers  to  simulate  a  damage  within  the  sample  (see  [3,  8]  for  graphical  representations). 
The  “damaged”  piece  of  aluminum  is  moved  from  one  layer  to  another  to  simulate  damages  within  the 
sample  at  different  depths,  and  the  length  of  the  damage  is  varied  by  producing  “gaps”  of  varying  size 
from  the  aluminum  plate  (the  thickness  of  the  damage  is  always  fixed  at  1mm).  As  a  means  of  inducing 
current  within  the  sample,  a  thin  sheet  of  copper  carrying  a  uniform  current  of  3  Amps  is  placed  above 
the  sample  on  top  of  a  thin  sheet  of  paper  (to  avoid  direct  physical  contact  between  the  sample  and  the 
conducting  sheet).  The  GMR  sensor  measures  the  amplitude  and  phase  of  the  magnetic  flux  density  across 
a  2in  line  (along  the  length  of  the  sample)  every  0.635mm.  The  data  is  then  filtered  through  a  lock-in 
amplifier  and  saved  to  a  file. 

The  experimental  data  did  not  resemble  Ansoft  simulations  from  our  model;  therefore,  we  chose  to 
snapshot  on  the  data  itself  in  forming  the  POD  basis  elements.  Furthermore,  it  was  necessary  to  filter  out 
the  background  noise  (data  for  a  sample  containing  no  damage)  in  order  to  obtain  a  representative  pattern 
in  the  data  as  a  function  of  the  damage  within  the  sample  (see  [3]  for  full  details).  In  filtering  out  the 
background  noise,  we  magnified  boundary  effects  from  the  edges  of  the  sample.  Therefore,  in  the  cost 
criterion,  we  only  consider  data  across  the  center  of  the  sample.  In  other  words,  we  use  the  cost  criterion 
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given  by 


-  5] 


(19) 


where  q  is  the  vector  containing  the  parameters  we  wish  to  estimate,  (q)  is  the  POD  approximation 
formed  using  snapshots  on  the  data  itself,  B2  1®  OMR  data  at  grid  points  Xj,  j  =  1, n  with  n  total  grid 
points  and  a  and  b  indicate  the  indices  of  the  grid  points  we  consider  in  our  cost  criterion. 

In  estimating  a  depth  of  d*  =  3mm  with  corresponding  length  T  =  1.5cm,  we  obtained  an  estimate 
of  d  =  2.9759mm  and  I  =  1.4612cm,  resulting  in  a  relative  error  for  depth  of  Ra  =  0.80%  and  length 
Ri  =  2.59%.  Again,  even  with  experimental  data,  the  results  were  accurate. 


4  Conclusion 

In  this  paper  we  presented  a  reduced  order  computational  algorithm  which  contributes  to  the  overall  field  of 
nondestructive  evaluation.  In  detecting  subsurface  damages,  there  is  a  need  for  a  fast  and  efficient  inverse 
problem  methodology.  The  reduced  order  POD  method  is  an  attractive  method,  because  it  allowed  us  to 
create  a  reduced  basis  of  less  than  10  basis  elements  in  the  trials  presented  while  still  capturing  99%  of  the 
total  energy  of  the  system.  This  results  in  an  accurate,  as  well  as  fast,  forward  algorithm. 

We  then  gave  an  overview  of  the  POD  methodology  in  the  context  of  subsurface  damage  detection 
or  parameter  identification.  We  discussed  the  implementation  of  the  inverse  problem  for  both  simulated 
results  as  well  as  experimental  results,  providing  examples  for  each  case.  Based  upon  the  results  presented, 
we  showed  that  using  either  simulated  data  or  experimental  data  on  either  a  one-parameter  estimation 
problem  or  a  two-parameter  estimation  problem,  we  achieved  accurate  results. 

However,  the  most  significant  finding  is  in  regard  to  reduction  in  computational  time  which  can  be  sum¬ 
marized  as  follows.  If  one  were  to  use  a  software  package  such  as  Ansoft’s  Maxwell  2D  Field  Simulator  to 
calculate  the  forward  problem  each  time  it  is  required  in  the  inverse  problem,  it  would  take  approximately 
5-10  minutes  for  a  single  forward  solve.  The  forward  solve  is  typically  required  anywhere  from  as  few 
as  10  times  to  as  many  as  500-1000  times.  Assuming  the  forward  algorithm  is  called  approximately  100 
times  with  an  average  of  7  minutes  for  each  forward  run,  the  inverse  problem  would  take  approximately 
42,000  seconds.  In  the  examples  provided,  the  entire  inverse  problem  ran  in  approximately  8-10  seconds. 
Consequently,  on  average  the  computational  time  is  reduced  by  a  factor  of  10^.  This  suggests  that  an  ap¬ 
propriate  sensing  device,  when  coupled  with  reduced  order  modeling  in  the  inverse  problem,  might  prove 
feasible  in  practical  damage  detection  applications. 
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